Cod diet: spatial analysis of cod and flounder stomach content in relation to density

Intro

Spatial analysis of cod and flounder stomach content in relation to density using sdmTMB

# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 12))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)

# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/stomach_fullness_spatial_cache/html")

For maps

# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

ggplot(swe_coast_proj) + geom_sf() 

# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'black', margin = margin()),
      strip.background = element_rect(fill = "grey90")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'black', margin = margin()),
        strip.background = element_rect(fill = "grey90")
      )
}

# Make default base map plot
plot_map_raster <- 
ggplot(swe_coast_proj) + 
  geom_sf(size = 0.3) +
  labs(x = "Longitude", y = "Latitude") +
  theme_facet_map(base_size = 14)

# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
  xy <- data.frame(ID = 1:length(x), X = x, Y = y)
  coordinates(xy) <- c("X", "Y")
  proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
  res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
  return(as.data.frame(res))
}

theme_clean <- function() {
  theme_minimal(base_family = "Barlow Semi Condensed") +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
          axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
          strip.text = element_text(family = "BarlowSemiCondensed-Bold",
                                    size = rel(1), hjust = 0),
          strip.background = element_rect(fill = "grey80", color = NA))
}

Read and plot data

# These data are for total and prey specific stomach models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   time_period = col_character(),
#>   quarter_fact = col_character(),
#>   pred_weight_source = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.

cod <- cod %>%
  mutate(year = as.integer(year),
         quarter = as.factor(quarter),
         depth2_sc = depth - mean(depth),
         saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  filter(year > 2014) %>%
  filter(!quarter == 2) %>% 
  drop_na(predfle_density_sc, predcod_density_sc) %>% 
  droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#>         converted 'quarter' from double to factor (0 new NA)
#>         new variable 'depth2_sc' (double) with 102 unique values and 0% NA
#>         new variable 'saduria_entomon_per_mass' (double) with 1,666 unique values and 0% NA
#>         new variable 'tot_prey_biom_per_mass' (double) with 7,108 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 102 unique values and 0% NA
#> filter: removed 4,936 rows (58%), 3,593 rows remaining
#> filter: no rows removed
#> drop_na: no rows removed

fle <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   pred_weight_source = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.

fle <- fle %>%
  mutate(year = as.integer(year),
         quarter = as.factor(quarter),
         depth2_sc = depth - mean(depth),
         saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  filter(!quarter == 2) %>% 
  drop_na(predfle_density_sc, predcod_density_sc) %>% 
  droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#>         converted 'quarter' from double to factor (0 new NA)
#>         new variable 'depth2_sc' (double) with 80 unique values and 0% NA
#>         new variable 'saduria_entomon_per_mass' (double) with 596 unique values and 0% NA
#>         new variable 'tot_prey_biom_per_mass' (double) with 1,699 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 80 unique values and 0% NA
#> filter: no rows removed
#> drop_na: no rows removed

# Plot data (with smoother)
ggplot(fle, aes(predfle_density_sc, tot_prey_biom_per_mass)) +
    geom_point() +
    stat_smooth(method = "gam", formula = y ~ s(x, k = 3))

  
ggplot(fle, aes(predcod_density_sc, tot_prey_biom_per_mass)) +
  geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3))

ggplot(cod, aes(predfle_density_sc, tot_prey_biom_per_mass)) +
    geom_point() +
    stat_smooth(method = "gam", formula = y ~ s(x, k = 3))

  
ggplot(cod, aes(predcod_density_sc, tot_prey_biom_per_mass)) +
  geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3))

Read the prediction grids

# And now read in pred_grid2 which has oxygen values at location and time and depth:
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_interactions/main/data/pred_grid2.csv")

# Clean
pred_grid2 <- pred_grid2 %>%
  mutate(year = as.integer(year)) %>% 
  drop_na(depth)

# Add ices_rect
pred_grid2$ices_rect <- ices.rect2(pred_grid2$lon, pred_grid2$lat) 

Fit sdmTMB models with densities as covariates to stomach content data. First make mesh

# Cod 
pcod_spde <- make_mesh(cod, c("X", "Y"), n_knots = 150, type = "kmeans", seed = 42)
plot(pcod_spde)

pfle_spde <- make_mesh(fle, c("X", "Y"), n_knots = 70, type = "kmeans", seed = 42)
plot(pfle_spde)

Fit models of total content

# Model total prey biomass 
# Cod 
mcod <- sdmTMB(
  data = cod, 
  formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
  time = "year", spatiotemporal = "IID", spatial = "on", mesh = pcod_spde, family = tweedie())

tidy(mcod)
#>                  term    estimate  std.error
#> 1            quarter1 -4.67446730 0.18369183
#> 2            quarter4 -4.62368659 0.15949784
#> 3 as.factor(year)2016  0.16511474 0.19693266
#> 4 as.factor(year)2017  0.30589646 0.20053012
#> 5 as.factor(year)2018  0.40552288 0.21650533
#> 6 as.factor(year)2019  0.27353933 0.22597431
#> 7 as.factor(year)2020  0.33937416 0.23827154
#> 8  predfle_density_sc -0.08755258 0.05474359
#> 9  predcod_density_sc  0.03191979 0.05142004
summary(mcod)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc
#> Time column: "year"
#> Mesh: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -4.67    0.18
#> quarter4               -4.62    0.16
#> as.factor(year)2016     0.17    0.20
#> as.factor(year)2017     0.31    0.20
#> as.factor(year)2018     0.41    0.22
#> as.factor(year)2019     0.27    0.23
#> as.factor(year)2020     0.34    0.24
#> predfle_density_sc     -0.09    0.05
#> predcod_density_sc      0.03    0.05
#> 
#> Dispersion parameter: 0.53
#> Tweedie p: 1.74
#> Matern range: 15.22
#> Spatial SD: 0.12
#> Spatiotemporal SD: 0.61
#> ML criterion at convergence: -9904.608
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

cod$resids_mcod <- residuals(mcod) # randomized quantile residuals
qqnorm(cod$resids_mcod); abline(a = 0, b = 1)

# Flounder
mfle <- sdmTMB(
  data = fle, 
  formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
  time = "year", spatiotemporal = "IID", spatial = "on", mesh = pfle_spde, family = tweedie())

tidy(mfle)
#>                  term     estimate  std.error
#> 1            quarter1 -5.054968610 0.27046794
#> 2            quarter4 -4.560112128 0.24061922
#> 3 as.factor(year)2016 -0.028909687 0.27451881
#> 4 as.factor(year)2017  0.301285213 0.26933035
#> 5 as.factor(year)2018  0.672111210 0.33586941
#> 6 as.factor(year)2019  0.620685733 0.28679830
#> 7 as.factor(year)2020  1.026916621 0.28115647
#> 8  predfle_density_sc  0.002925394 0.09706364
#> 9  predcod_density_sc  0.055584176 0.11765389
summary(mfle)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc
#> Time column: "year"
#> Mesh: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -5.05    0.27
#> quarter4               -4.56    0.24
#> as.factor(year)2016    -0.03    0.27
#> as.factor(year)2017     0.30    0.27
#> as.factor(year)2018     0.67    0.34
#> as.factor(year)2019     0.62    0.29
#> as.factor(year)2020     1.03    0.28
#> predfle_density_sc      0.00    0.10
#> predcod_density_sc      0.06    0.12
#> 
#> Dispersion parameter: 0.12
#> Tweedie p: 1.50
#> Matern range: 39.88
#> Spatial SD: 0.57
#> Spatiotemporal SD: 0.54
#> ML criterion at convergence: -4667.908
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

fle$resids_mfle <- residuals(mfle) # randomized quantile residuals
qqnorm(fle$resids_mfle); abline(a = 0, b = 1)

Fit models of saduria content

# Model Saduria contents 
# Cod
mcodsad <- sdmTMB(
  data = cod, 
  formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
  time = "year", spatiotemporal = "IID", spatial = "on", mesh = pcod_spde, family = tweedie())

tidy(mcodsad)
#>                  term    estimate std.error
#> 1            quarter1 -8.87619104 0.7663188
#> 2            quarter4 -9.08924932 0.7214733
#> 3 as.factor(year)2016 -0.82029185 0.7820977
#> 4 as.factor(year)2017 -1.36653144 0.7993448
#> 5 as.factor(year)2018 -0.26043020 0.8378119
#> 6 as.factor(year)2019 -0.88611599 0.9053763
#> 7 as.factor(year)2020 -0.01723844 0.8681083
#> 8  predfle_density_sc -0.60227548 0.3012112
#> 9  predcod_density_sc  0.39255981 0.2620735
summary(mcodsad)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc
#> Time column: "year"
#> Mesh: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -8.88    0.77
#> quarter4               -9.09    0.72
#> as.factor(year)2016    -0.82    0.78
#> as.factor(year)2017    -1.37    0.80
#> as.factor(year)2018    -0.26    0.84
#> as.factor(year)2019    -0.89    0.91
#> as.factor(year)2020    -0.02    0.87
#> predfle_density_sc     -0.60    0.30
#> predcod_density_sc      0.39    0.26
#> 
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> Matern range: 36.25
#> Spatial SD: 1.88
#> Spatiotemporal SD: 1.56
#> ML criterion at convergence: -1024.766
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

cod$resids_mcodsad <- residuals(mcodsad) # randomized quantile residuals
qqnorm(cod$resids_mcodsad); abline(a = 0, b = 1)

# Flounder
mflesad <- sdmTMB(
  data = fle, 
  formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc,
  time = "year", spatiotemporal = "IID", spatial = "on", mesh = pfle_spde, family = tweedie())

tidy(mflesad)
#>                  term   estimate std.error
#> 1            quarter1 -8.7008421 1.0345886
#> 2            quarter4 -7.9176420 0.9779841
#> 3 as.factor(year)2016  0.7217416 0.7990698
#> 4 as.factor(year)2017  0.5084692 0.7999960
#> 5 as.factor(year)2018  1.8324164 0.9723237
#> 6 as.factor(year)2019  0.1362148 0.8635915
#> 7 as.factor(year)2020  1.7146474 0.8095086
#> 8  predfle_density_sc -0.1720328 0.3376522
#> 9  predcod_density_sc  0.6932675 0.4896959
summary(mflesad)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc
#> Time column: "year"
#> Mesh: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -8.70    1.03
#> quarter4               -7.92    0.98
#> as.factor(year)2016     0.72    0.80
#> as.factor(year)2017     0.51    0.80
#> as.factor(year)2018     1.83    0.97
#> as.factor(year)2019     0.14    0.86
#> as.factor(year)2020     1.71    0.81
#> predfle_density_sc     -0.17    0.34
#> predcod_density_sc      0.69    0.49
#> 
#> Dispersion parameter: 0.14
#> Tweedie p: 1.48
#> Matern range: 66.65
#> Spatial SD: 2.44
#> Spatiotemporal SD: 1.18
#> ML criterion at convergence: -1440.375
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

fle$resids_mflesad <- residuals(mflesad) # randomized quantile residuals
qqnorm(fle$resids_mflesad); abline(a = 0, b = 1)

Plot effect sizes

mcod_ef <- tidy(mcod, effects = "fixed", conf.int = TRUE) %>%
  filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>% 
  mutate(species = "Cod", response = "Total prey")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

mfle_ef <- tidy(mfle, effects = "fixed", conf.int = TRUE) %>%
  filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>% 
  mutate(species = "Flounder", response = "Total prey")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

mcodsad_ef <- tidy(mcodsad, effects = "fixed", conf.int = TRUE) %>%
  filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>% 
  mutate(species = "Cod", response = "Saduria")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

mflesad_ef <- tidy(mflesad, effects = "fixed", conf.int = TRUE) %>%
  filter(term %in% c("predfle_density_sc", "predcod_density_sc")) %>% 
  mutate(species = "Flounder", response = "Saduria")
#> filter: removed 7 rows (78%), 2 rows remaining
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'response' (character) with one unique value and 0% NA

plot_df <- bind_rows(mcod_ef, mfle_ef, mcodsad_ef, mflesad_ef) %>% 
  mutate(term = ifelse(term == "predcod_density_sc", "Cod density", "Flounder density"))
#> mutate: changed 8 values (100%) of 'term' (0 new NA)

# Plot effects
plot_df %>% 
  ggplot(., aes(term, estimate, color = factor(species))) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.75) +
  geom_point(size = 3, position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.2), size = 1) +
  scale_color_brewer(palette = "Dark2", name = "Species") +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme_light(base_size = 12) + 
  facet_wrap(~response) +
  theme(legend.position = "bottom") +
  NULL 

ggsave("figures/stomach_content_effect_size.png", width = 6.5, height = 6.5, dpi = 600)

Calculate and plot marginal effects

Cod

# Flounder density
nd_cod_fle <- data.frame(predfle_density_sc = seq(min(cod$predfle_density_sc),
                                                  max(cod$predfle_density_sc), length.out = 100))

nd_cod_fle$year <- 2018L
nd_cod_fle$predcod_density_sc <- 0
nd_cod_fle$depth_sc <- 0
nd_cod_fle$quarter <- factor(4)

sad_pred_cod_fle <- predict(mcodsad, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
tot_pred_cod_fle <- predict(mcod, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)

# Cod density
nd_cod_cod <- data.frame(predcod_density_sc = seq(min(cod$predcod_density_sc),
                                                  max(cod$predcod_density_sc), length.out = 100))

nd_cod_cod$year <- 2018L
nd_cod_cod$predfle_density_sc <- 0
nd_cod_cod$depth_sc <- 0
nd_cod_cod$quarter <- factor(4)

sad_pred_cod_cod <- predict(mcodsad, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
tot_pred_cod_cod <- predict(mcod, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)

Flounder

# Flounder density
nd_fle_fle <- data.frame(predfle_density_sc = seq(min(fle$predfle_density_sc),
                                                  max(fle$predfle_density_sc), length.out = 100))

nd_fle_fle$year <- 2018L
nd_fle_fle$predcod_density_sc <- 0
nd_fle_fle$depth_sc <- 0
nd_fle_fle$quarter <- factor(4)

sad_pred_fle_fle <- predict(mflesad, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
tot_pred_fle_fle <- predict(mfle, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)

# Cod density
nd_fle_cod <- data.frame(predcod_density_sc = seq(min(fle$predcod_density_sc),
                                                  max(fle$predcod_density_sc), length.out = 100))

nd_fle_cod$year <- 2018L
nd_fle_cod$predfle_density_sc <- 0
nd_fle_cod$depth_sc <- 0
nd_fle_cod$quarter <- factor(4)

sad_pred_fle_cod <- predict(mflesad, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
tot_pred_fle_cod <- predict(mfle, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)

Plot marginal effects

# Saduria models
p1 <- ggplot(sad_pred_fle_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  labs(y = "Saduria in flounder stomach [g/g]", x = "")

p2 <- ggplot(sad_pred_fle_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  labs(y = "", x = "")

p3 <- ggplot(sad_pred_cod_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  labs(y = "Saduria in cod stomach [g/g]", x = "Scaled flounder density")

p4 <- ggplot(sad_pred_cod_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  labs(y = "", x = "Scaled cod density")

p1 + p2 + p3 + p4 + plot_layout(ncol = 2)

ggsave("figures/marginal_effects_saduria.png", width = 6.5, height = 6.5, dpi = 600)

# Total stomach content models
p5 <- ggplot(tot_pred_fle_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  labs(y = "Prey biomass in flounder stomach [g/g]", x = "")

p6 <- ggplot(tot_pred_fle_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  labs(y = "", x = "")

p7 <- ggplot(tot_pred_cod_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  labs(y = "Prey biomass in cod stomach [g/g]", x = "Scaled flounder density")

p8 <- ggplot(tot_pred_cod_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  labs(y = "", x = "Scaled cod density")

p5 + p6 + p7 + p8 + plot_layout(ncol = 2)

ggsave("figures/marginal_effects_total.png", width = 6.5, height = 6.5, dpi = 600)

knitr::knit_exit()